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Abstract 

Kinetic energy driven phase transitions in Bose superfluids occur at low values of the 
repulsion when the values of the next-to-nearest and next-to-next-to-nearest hopping term 
attain certain critical values, resulting in alterations in the wave vector of the condensate. 
We map out the space of possible phases allowed by particular forms of the single-particle 
energy dispersion in the superfluid state, noting the appearance of a new phase, and 
examine in more detail the effects of additional repulsive terms on the form of the con- 
densate wavefunction. We also examine the effect of these additional hopping terms on 
the formation of inhomogeneities in the condensate. 

PACS numbers: 
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I. INTRODUCTION 



Ever since the phase transition (QPT) from a superfluid hquid to a Mott in- 
sulator was observed in atomic gases trapped in an optical lattice [1], there has 
been considerable experimental and theoretical activity geared towards the explo- 
ration of the physics of this transition as described by the Bose-Hubbard model 
This phenomenon has also long been of interest in the study of phase transitions 
of charged Bose liquids as observed in the context of granular superconductors and 
Josephson junction arrays [3], in addition to preformed real-space electron pairs 
such as bipolarons 

The superfluid-Mott transition occurs only at fixed densities in a homogeneous 
system when the number of bosons is commensurate with the number of lattice 
sites. At intermediate densities a transition may occur to a density wave superfiuid 
jsl, Q] (similar to the putative supersolid phase of ^He \}\) which has a tendency 
towards phase separation 0, 1^, Q. The transitions require that both the repulsive 
interaction and the number of bosons is sufficiently large: the critical value of 
this repulsion for commensurate density has 3een determined numerically for two- 



dimensional square lattices 
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Alexandrov and Thomas [16] have shown that the inclusion of modifications 



to the kinetic energy term of the tight-binding Gross-Pitaevskii |l7l . Il8l | equation 
(a mean field description of the Bose-Hubbard Model [19]), along with additional 
longer-range repulsive terms gives rise to condensate stripes or checkerboards (as 
anticipated for d-wave superconductors in {2^) with wave- vectors that are not nec- 
essarily commensurate with the lattice period. Unusually, it is the values of the 
additional hopping parameters that drive these 'kinetic' phase transitions between 
the homogeneous phase and the striped/checkerboard phases, and not the repul- 
sive potentials - though the repulsion does play a role regarding whether stripes 
or checkerboards are formed. In the present paper, we expand on those results, 
examining these phenomena in more detail. 

At first glance, it might not seem that one could expect a phase transition in 
the ground state of a near-ideal Bose gas whose inter-particle repulsion is small 
compared to their kinetic energy. In particular, this absence of any significant 
effects is suggested by the rather strong constraints placed on the form of the 
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30Sonic ground state wave function by the requirement that it posses no nodes 
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221], which imphes that the values of the hopping terms for the tight-binding 



form of the s-wave ground state must be such that they approximate a kinetic 
dispersion with a single minimum located at the zero wave-vector. There are a 
number of subtleties, however, which appear to indicate that this is not all that can 
be said on the matter. 

Firstly, there exist some cases (such as that of bipolarons) where the condensation 
occurs in the p-band (or the d-band in the case of single polarons), and thus one's 
choice of hopping parameters is no longer limited by the no-node restriction, since 
this is not the ground state of the system as a whole. Another potential example 
of this is the incommensurate, met ast able superfluidity of atoms in the flrst excited 
)and of a double- well optical lattice that has been predicted by Stojanovic et al 



23[ |. since it has been shown that the lifetimes of excited atoms in p-wave orbitals 



in optical lattices may be extended by an unexpectedly large amount [2J, |25 |. 

Secondly, the assumptions of this no-go theorem [22] are violated in the presence 
of signiflcant many-body eflFects. This can be explicit, as in the case of the non- 
linear potential in the GPE (which renders the Hamiltonian non-self- adjoint), or 
implicit and therefore hidden in the form of the hopping terms such that they no 
longer approximate the continuum Laplacian operator of the kinetic term. This last 
can occur, for example, as a result of interactions with other systems that produces 



long-range dispersive effects such as in Keverekidis et a/.'s [26|] study of the dynamics 
of discrete breathers in a coupled linear and non-linear system, or such as Feshbach 
resonance-type effects that result in the p-wave state becoming the system ground 



state 



271 . 1281] . Similarly, Larson and Martikainen have shown that the interaction 



between internal and external degrees of freedom in an ultra-cold gas of atoms can 
give rise to a transition to a ground state with nodes as the detuning of the atom 



fleld is adjusted 



291 ]: when written in terms of an effective Bose-Hubbard model. 



this new ground state corresponds to a negative value of the nearest-neighbour 
hopping parameter jsol - that is, the dispersion has a p-wave form, characterised 
by an 'antiferromagnetic' ordering of the condensate wave function. 

Given the above considerations, it is worth exploring the behaviour of bosonic 
superfluids on a lattice in the GPE limit for all relevant values of the additional 
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hopping parameters, despite possible difficulties in the strai 



replication of these values on certain optical lattices [3j|. Similar analytic and 



ghtforward experimental 



numerical analyses have been performed in the case of the full Bose Hubbard model 



331], though not for precisely identical parameters. In what follows, we outline 



our formalism and numerical method (^III), expand on the notion of kinetic energy 
driven phase transitions ( ^TTTl) first described in {igI and noting the appearance of 
a new phase not described earlier, before discussing our numerical results in the 
context of this theory (^TVland fVT). 



II. THE TIGHT-BINDING GROSS-PITAEVSKII EQUATION 

We begin with a system whose ground state is described by a Gross-Pitaevskii- 
type (GP) equation [l7|, ll8(] including lattice, V{r)^ and interaction, t/(r), poten- 
tials, 



-^ + Vir)-^i + j dv'U{v-v'mv')\^ m = 0, (1) 

where the condensate wave- functions ip{r) play the role of the order parame- 
ter m is the boson mass, and /x is the chemical potential. '0(r) is normalised 
through J rfr|'0(r)p = A^^, where N^, is the number of bosons determined by the 
value of /X. 

To solve equation ([1]) variationally on a lattice with a large number of sites, 
we make use of a complete set of orthogonal Wannier (site) functions w{r). Writing 
the order parameter in terms of these as ip{r) = "^^m^^^i^ ~ reduce 
([T]) to a discrete set of equations for the amplitudes at each lattice site, 

- 5][t(m - m') + ^S^,n.']<Pm + Yl ^mn"'C'0n0m = 0. (2) 
m m,n,n^ 

Here t(m) = J (iri(;*(r)[— /i^V^/2m + V{r)]w{r — m) is the hopping integral, and 
^mn"' ^11 drdr^U{r — r^)t(;*(r — m^)t(;*(r^ — n^)w{r' — n)w{r — m) is the matrix 
element of the interaction potential. 

In the liquid regime far away from the superfiuid-Mott transition, the hopping 
integrals dominate over the interaction matrix elements, so that we need only keep 
the density-density interactions, U^^^' ^ U^^S^ ^n' Sn,n' • Solutions to the GP equa- 
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tion ([2D are then equivalent to the minima of the energy functional, 

1 



m,n 



t(m - n) + flS^^n - o^mnrmV^n 



2" 



(3) 



' mn 



Rescaling the order parameter as = n^^^fm and the interaction as 
Unm/'f^ yields a universal functional, E{f^) = E{(f)^)/n that is independent of 
n = Ni)/N ^ the number of bosons per lattice valley: 



E{u = - E 



m,n 



t(m - n) + n5ra,n - ^UmnfLfn 



f:fm. (4) 



Importantly, this shows that solutions for different particle densities can be mapped 
on to each other by a simple rescaling of the interaction. 




FIG. 1: Diagram showing the different hopping parameters from a given site in the forward 
right direction on a 2D lattice. 

In what follows we examine the effects of the variation of the values of t2 and ^3, 
the next-to-nearest neighbour and the next-to-next-to-nearest neighbour parame- 
ters respectively (Fig. ([T])), for various small values of the onsite (t^oo), nearest 
neighbour (uqi) and next-to-nearest neighbour (1^02) repulsive potentials. 

Numerical results were obtained through the following procedure: we generated 
random starting values of the amplitudes /m (where /m ^ 3? because there is no 
current) on a lattice with sites and periodic boundary conditions. The functional 
dH) was minimised for a given value of fi using the NAG conjugate gradient opti- 
misation routine E04DGF. When calculated as part of this procedure, the kinetic 
portion of the functional was calculated in momentum space following a Fourier 
transformation of the amplitudes at each lattice site. This was then inverse Fourier 
transformed back into Wannier space where the non-linear potential portion of the 
functional was then applied. Note that the algorithm tends to fail in the absence of 
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a small on-site repulsion. We fine-tuned the value of /x using the Van Wijngaarden- 
Dekker-Brent algorithm [34^ so that the constraint l/mP = ^ was imposed; 
typically /x is determined to an accuracy of 10"'' — 10~^, corresponding to a relative 
error of order 10~^ — 10"'' in the enforcement of the constraint, which is sufficient 
for our purposes. Unless otherwise stated, a square 2D lattice of = 50^ has been 
used for the plots of the condensate function; simulation on a variety of lattice sizes 
seems to indicate that this is most likely large enough to avoid significant finite 
size effects. However, 25 x 25 lattices seem to give good qualitative indications of 
the phase at a given point, which should be largely independent of the lattice size. 
Due to the relatively quick speed with which the minimisation algorithm runs for 
lattices of this size, the results in Figure [I2] pertaining to the scanning of the t2 
and ts parameter spaces were generated on these. Some cross-checks for certain 
parameter sets were performed on 50 x 50 lattices; when in doubt, these were used 
to determine the phase of the condensate. 

It is worth noting that one sometimes finds in the ordered phases that a pair 
of parallel 'kinks' (or a quartet, if they are ordered diagonally) has formed in the 
condensate, and that these are often slightly more stable than the ordered state 
alone. We suspect that these are artefacts of the periodic boundary conditions 
and the finite lattice size, perhaps made frequent by our use of a random start 
instead of a trial wave-function; firstly because at small repulsions, the 'kinks' 
vanish - or become much less common ^ when the lattice size is set to 50 x 50, 
and secondly because they always appear in pairs or quartets so as to satisfy the 
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boundary conditions. We therefore neglect them; even if they are not artefacts the 
presence of 'kinks' in the order of the condensate are not as important to us as the 
nature of that order, which is the focus of our study, and they are quahtatively 
different from the meandering domain waUs discussed in §IV[ 

In addition, on the 25 x 25 lattices one finds that in the commensurate phases 
there is a dislocation; this is an artefact both of the periodic boundary conditions 
and there being an odd number of lattice sites in each direction, since this entails 
that the wavefunction at one end is out of phase with the wavefunction at the other 
end. This is not an obstacle to the correct diagnosis of the phase of the system, 
but for this reason any figures presented displaying numerical calculations of the 
ground state wavefunction have been performed using 50 x 50 lattices. 




FIG. 3: Dispersion relations for a) t2 = —0.5 and b) = —0.25 in the conventional case, 
and c) t2 = —0.5 and d) ts = —0.25 in the unconventional case . 
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III. KINETIC ENERGY DRIVEN PHASE TRANSITIONS 



Provided that the repulsion is weak, the ground state of the system is determined 
by the form of the kinetic contribution to the Hamiltonian, which in the tight 
binding hmit incorporates the effects of the periodic lattice into the values of the 
hopping parameters t\ 

E(k) = 4to — 2ti{co8{kx) + cos(/c^)) — 4^2 co8{kx) cos(/c^) — 2ts{co8{2kx) + cos(2/c^)), 

(5) 

where for the purposes of our present discussion we take to = (its precise value 
is of little relevance here). The basic form of the band is given by the sign of ti 
(whose absolute value is taken to be unity here and thereafter): an conventional 
dispersion has ti positive and a single minimum at the F point and a unconventional 
dispersion has ti negative and four minima located in the corners of the Brillouin 
zone (Figure [2]). At small values of the repulsion the condensate wavefunction in 
the former case is flat, in the latter case it has a staggered (alternating) sign with 
each lattice site. In both cases the condensate density is homogeneous. 




FIG. 4: Dispersion relations for the unconventional case with a) t2 = —0.8 (C phase) and 
b) ts = -0.25 (D phase). 

We restrict our discussion to cases where t2^ts < 0, since this appears to be where 
many of the phenomena of interest reside. As discussed in [16], we may modify the 
behaviour of the system by decreasing the values of the hopping parameters t2 and 
t^. Varying t2 only or only we discover that the value of the inverse effective mass 
is decreased until at certain critical values ^ —0.5 and ^ —0.25) it becomes 
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zero (this corresponds to a flattening of the dispersion in the vicinity of the minima 
- see Figure [3]). As our discussion in ^indicates, in the conventional case we should 
take these to be the lower limit of the allowed values of those hopping parameters, 
since beneath these values the ground state of the system develops nodes, and 
this is unphysical for the ground state of a bosonic system. In the unconventional 
case, however, we take this to be an indication of a phase transition where the wave- 
vect or /inverse period of the condensate modulation is the order parameter; here the 
nodes of the system change location. When ^2 < the minima are now located 
at at ki = (tt, 0), k2 = (— tt, 0), ks = (0, tt) and k4 = (0, — tt); this is a flrst-order 
phase transition due to the discontinuous change in the location of the minima. 
We shall call this the Commensurate phase (C phase) due to the appearance 
of modulations in the ground state that are commensurate with the lattice period. 
When ts < t^c^ these are located at ki = /c(l,l),k2 = 1,1), ks = 1,— 1) 
and k4 = /c(l, — 1), where is a constant factor determined by the value of that 
is equal to tt (unconventional case) or (conventional) at and tends towards 
7r/2 as is decreased |l6|. This results in modulations in the ground state that 
are aligned diagonally but that are not commensurate with the lattice period: we 
therefore denote this as the Diagonal phase (D phase). Due to the continuous 
change in the order parameter, this is a second-order phase transition. Examples 
of dispersion relations within these phases may be found in Figure [H We call the 
phase transitions 'kinetic energy driven' as it is the form of the tight-binding kinetic 
dispersion relation that drives them. 

The situation when both t2 and ts are varied together is somewhat more complex, 
as varying both has an effect on the location of the transitions to the various phases. 
In addition, we may see that there exists an additional phase. In the conventional 
case, we have four minima located at ki = (/,0), k2 = — (/,0), ka = (0,/) and 
k4 = 7r(0, /) where / is dependent on the values of the hopping parameters and 
tends towards tt if t2 alone is increased. The unconventional case is equivalent to a 
displacement of the centre of the conventional Brillouin zone by (7r/2,7r/2) and so 
there are eight minima instead of four, located at k^ = 7r(l, ±/), k2 = — 7r(l, ±/), 
kg = 7r(±/, 1) and = — 7r(±/, 1), where / is again dependent on the values of the 
hopping parameters and tends towards if ^2 alone is increased. Figure [5] shows the 
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-1.9 




FIG. 5: Low-energy portions of dispersions in the unconventional case with = —0.2 
and a) t2 = —0.4, b) = —0.5, c) t2 = —0.8 and d) t2 = —0.9, showing the rotation of 
the minima. 



evolution of the wave vectors for ts = —0.2 as t2 is varied - in this case, the wave- 
vectors are rotated away from the diagonal axis of the D phase, until they finally 
reach the vertical and horizontal axes of the C phase. We therefore designate this 
to be the Rotated phase (R phase). 

Examination of the form of the dispersion relation for different values of t2 and 
ts allows us to classify various regions of the parameter space according to these 
phases. This is done in Figure [6l 

One method of realising this behaviour in optical lattices is through the use 
of shallow lattice potentials, such that the localisation of the Wannier functions 
describing the lattice is poor. Stojanovic et al {2^ predict incommensurate con- 
densate wave vectors for atoms in the meta-stable p-band of a double- well lattice for 
sufficiently small values of the lattice potential; the example given of a correspond- 
ing dispersion indicates that this is within the region of our D-phase. Note that in 
most cases involving optical lattices, conventional type dispersions with more than 
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FIG. 6: The form of the dispersion on the t2 — plane. White signifies the ordinary 
homogeneous phase, gray the D phase, black the C-phase and gold the R-phase. A 
divided circle indicates a dispersion intermediate between two phases, as in Figure [5] a). 

one zero are likely unphysical if we are in the s-wave ground state; in this case, 
Figure [6] displays possible bounds on the allowed values of t2 and ^3. 

In the next sections, we shall discuss the effects of the presence of different forms 
of potential on the condensate wave function for different values of the hopping 
parameters. 




FIG. 7: Inhomogeneity in the density for the unconventional case when i^oo = 4.0 and 
additional hopping terms are zero. 



11 



Simulation Type 


Approximate Uqq 


Gutzwiller-type Monte Carlo[12] 
Extended Gutzwiller-type Monte Carlofl3. 151 
Green's Function Monte Carlo[ll. 141 


23.32 
20.6 
16.7-17.0 



TABLE I: Numerical results showing the location of the superfluid-Mott insulator tran- 
sition for the Bose-Hubbard model at an average atomic density per site of n = 1 on a 
square, two-dimensional lattice. 



IV. INTERMEDIATE VALUES OF uqq AND ADDITIONAL HOPPING 
PARAMETERS 

What can be considered to be intermediate values of Uqq^ In Table [U we show var- 
ious critical values of Uqq resulting from numerical simulations of the Bose-Hubbard 
model on a 2D square lattice with commensurate density. In what follows, we take 
'intermediate' to mean any reasonably large {uqq > 1) value of the repulsion which 
is also below i^qq, and is therefore within the range of validity of the GPE. 

For sufficiently large (but still intermediate in the above sense) values of the 
onsite repulsion uqq^ we observe that the homogeneous ground-state becomes inho- 
mogeneous (Figure [71), with long, meandering domain walls. This is likely due to 
the repulsion becoming dominant over the periodic potential of the lattice, which 
in the tight-binding limit is incorporated into the kinetic portion of our Hamilto- 
nian. This effect in a bosonic system is reminiscient of a similar effect observed in 
fermionic systems (s^. 

An interesting feature of our model is the manner in which the critical value of 
the on-site repulsion at which this instability appears, i^qo' altered by the value 
of ^2- We have mapped out the approximate location of this transition through 
numerical simulation on 50 x 50 lattices for various values of Uqq and ^2, with ts 
and the other repulsions set to zero; our results are displayed in Figure El For 
the conventional (first discussed in jl6|) and unconventional cases, we see that u!^ 
decreases as ^2 decreases from towards —0.5, which correlates with the decrease 
in depth of the minimum of the dispersion relation. At t2 = —0.5, it is likely that 
t^QQ = 0. For the conventional case, decreasing t2 further is normally unphysical. 
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b) 

FIG. 8: Onset of inhomogeneity as uqq and t2 are varied. The gray regions are homoge- 
neous, the black inhomogeneous, and gold indicates the C phase, which has homogeneous 
density. Circles correspond to simulations whose resulting ground state is given by their 
color, a) displays the conventional case when restricted to hopping values giving a single 
zero (taken from fl6^), and b) displays the unconventional case. 



but for unconventional bosons this is not an issue; and for t2 < —0.5 we see that 
t^oo begins to increase once more until it returns to its original value at t2 = — 1- 

The behaviour of i^qq as t2 is varied corresponds closely with the height of the 
potential wall that separates one minimum from another (in the conventional case, 
that would be the height of the potential that separates the minimum at the F point 
of the first Brillouin zone from the minima at the F points of the neighbouring Bril- 
louin zones). As t2 is decreased, the height of this potential wall is also decreased, 
until it becomes fiat at t2 = —0.5, as in Figures [3] a) and b). As the height of the 
wall is decreased, the easier it is for the non-linear repulsion to overcome the kinetic 
energy. 
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E(k) 




FIG. 9: Plots of a) the dispersion and the condensate wavefunctions and densities respec- 
tively for b) and c) uqq = 0.02, d) and e) uqq = 0.02 and uqi = 0.01, f) and g) uqo = 0.02 
and uo2 = 0.02 in the unconventional case where ^2 = and ts = —0.4. Colours are 
assigned to contour plots b)-g) according to value of the far corner of a square. 
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FIG. 10: Plots of a) the dispersion and the condensate wavefunctions and densities re- 
spectively for b) and c) uqq = 0.02, d) and e) uqq = 0.02 and uqi = 0.01, f) and g) 
uqo = 0.02 and uo2 = 0.02 in the unconventional case where t2 = —0.6 and ts = —0. 
Colours are assigned to the contour plots for b)-g) according to value the point at the 
foremost corner of a square. No scale is given for c) and g) since the changes in the size 
of the density are <^ 1 and due to numeriq^ noise. 



For the unconventional case, further decreasing t2 results in a growing potential 
wall between the minima of the C-phase (compare Figure [3] c) with Figure [Ha), for 
example). As one would expect from the preceding argument, this entails that ul^ 
must increase. 

V. MODULATIONS FOR SMALL POTENTIALS 

We have performed the numerical minimisation procedure of ^for a number of 
values of the hopping parameters and small values of the repulsions. Figures M, flOl 
and [TT] shows results for parameter sets within the D, C and R phases respectively, 
where we have chosen a unconventional type dispersion since this is likely to be the 
most physically relevant. 

These results make the effects of the various forms of repulsion clear. In the 
non-homogeneous phases, a small repulsion causes the condensate wavefunction 
in Wannier space to take the form of a superposition respecting parity and time 
reversal symmetry, for which there are three possibilities, /m oc cos(ki • m), /m oc 
cos(k2 • m) and /m oc cos(ki • m) ± cos(k2 • m). The first two correspond to a 
'striped' condensate and the last corresponds to a 'checkerboard' condensate. Our 
figures show that a finite value of the nearest-neighbour repulsion uqi is needed for 
the selection of a checkerboard wavefunction in the C phases, whereas in the D phase 
a finite value of the next-nearest-neighbour repulsion U02 is required. Otherwise, 
we have incommensurate striped wavefunctions in the D phase, and commensurate 
striped wavefunctions resulting in a homogeneous condensate density in the C phase. 
The R-phase follows the same selection criteria for the form of its wave-function as 
does the C-phase, however using the k^2 instead. In addition, the wave-functions 
of both stripe and checkerboard phases are staggered due to the portions of the 
wave vectors that are located at the edges of the Brillouin zone; an conventional 
type dispersion would lack this staggering. 

It should be noted that the observed real-space modulations in our figures refiect 
the beating of the condensate period with that of the lattice. 

However, for combinations of values of the various repulsive terms, the relation- 
ship is not as straightforward. Figure [T2l shows scans of a region in the t2 — ts plane 
(carried out on 25 x 25 lattices with cross-checks on 50 x 50 lattices as discussed 
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FIG. 11: Plots of a) the dispersion and the condensate wavefunctions and densities re- 
spectively for b) and c) uqq = 0.02, d) and e) uqq = 0.02 and uqi = 0.01, f) and g) 
uqo = 0.02 and uo2 = 0.02 in the unconventional case where t2 = —0.6 and ts = —0.1. 
Colours are assigned to the contour plots for b)-g) according to value of the point at the 
foremost corner of a square. 
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FIG. 12: Scans of a region of the t2 — h plane with a) uqq = 0.2, uqi = 0.03, uo2 = 
0.02, h)uoo = 0.2, uoi = 0.01, uo2 = 0.02, c)uoo = 0.01, uqi = 0.03, uo2 = 0.02, d) 
uqo — 0.01, Uqi — 0.01, uo2 — 0.02, e) uqq — 0.2, ^xqi — 0.0, ixo2 — 0.02, and f) uqq — 
0.01, Uqi = 0.01, = 0.0. Filled white circles indicate a homogeneous condensate 
wave function, empty circles indicate a 'mixed' case, diagonal half-circles and crossed 
circles indicate diagonal stripes and checkerboards respectively, and vertical half-circles 
and crosses indicate stripes and checkerboards aligned with the axes of the lattice. In the 
latter cases, shading indicates that the stripe or checkerboard is incommensurate with 
the lattice. 



in ^TTl) for various values of the repulsive terms in order to obtain an exploratory 
overview of the behaviour of the system in the unconventional case. The potentials 
for these scans were selected with the following in mind: for t2 < —0.5 and ^3 = 0, 
we typically seem to require Uqq > Auqi in order to move from a checkerboard state 
to the commensurate striped state. 
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As a general rule, it seems that in order for a diagonal checkerboard to emerge, 
Uo2 > Uqo^ and that for an aligned (R and C phase) checkerboard state to emerge 
we must have Uqi > Uo2^ though this last condition does not seem to prevent the 
formation of a diagonal checkerboard in much (but not all) of the D phase. The 
> ^2 ^ —0.2, ts = —0.4 region appears to remain striped in the presence of a 
finite Uqi potential regardless of the value of Uo2- We may also observe the presence 
of 'mixed' ground states at certain values of ^2 and ^3; these correspond to the 
potential 'wall' separating diflFerent minima being shallow. This likely entails that 
in these areas the behaviour of the system is governed both by the potential terms 
and by the hopping terms in a complex fashion; this interplay would be worth 
further investigation. 

The tendency of the ^2 = —0.8, ts = —0.4 region to exhibit behaviours char- 
acteristic of both D and R phases is likely due to its intermediate character, this 
behaviour is also be apparent in the t2 = —0.4, = —0.2 region. However, the 
latter (when not in a 'mixed' phase) displays D-like behaviour when the former 
display R-like behaviour, and vice-versa. It is possible that there are slopes or 
other differences in the valleys of the dispersion of the their that are too small to 
be observed on the scale of our plots of the dispersion that may account for this. 

We are confident that our algorithm has in most cases located the appropriate 
ground states of the energy landscape (subject to the remarks in ^TT]), since they 
display behaviours consistent with what we would expect from Figure [6l We should 
note, however, that in the vicinity of the transition to a mixed phase the minimum 
corresponding to an ordered phase could be very close to that corresponding to 
a mixed phase, and the algorithm may have some difficulty in locating the true 
minimum. Therefore the results for ^2 = —0.2, = —0.2 should perhaps be treated 
with some caution. 

We might add that given the sensitivity of the form of the condensate wave- 
function to that of the potential, it would be interesting to observe the effects of 
the inclusion of a long range interaction - as in the case of a dipolar gas (such as 
in 0, 37, Q for example) or a magnetic ffeld in a charged Bose liquid 39, 40]. 
The latter would certainly be an interesting line of enquiry, given their role in a 
recently proposed explanation |4ll] of the anomalous magneto-oscillations recently 
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observed in cuprate superconductors [42|, |43[]. In addition, simulations of the full 



33] 



quantum Bose Hubbard model with finite positive nearest-neighbour hopping 
indicate possible quantum phase transitions to supersolid phases with finite values 
of the wave-vector. It will certainly prove worthwhile to determine any connec- 
tion between the transitions found here and those that might exist in the full Bose 
Hubbard model with negative values of the nearest-neighbour and next-to-nearest 
neighbour hopping parameters. 



VI. CONCLUSION 



In the foregoing, we have calculated numerical solutions of the Gross-Pitaevskii 
equation in the tight-binding limit with additional hopping parameters and non- 
linear repulsive terms. Building on the results of [16], we have described in more 
detail the mechanism of kinetic energy driven phase transitions in the superfiuid 
phase, locating a new phase in the process. We have also found that the value of Uqq 
at which one might expect an instability towards an inhomogeneous background to 
develop is altered by the variation of t2- This effect should be visible even in the 
case of conventional dispersions. Finally, we have also presented results showing the 
eflFects of the repulsive terms on the modulation of the condensate in the various 
phases, as well as performing an exploratory mapping of the parameter space for 
various values of the repulsive terms. 

We have discovered a new phase characterised by a wave vector which is incom- 
mensurate with the lattice in one direction and staggered in the other. In addition, 
we have found that the value of Uqq at which the condensate becomes unstable is 
affected by form of the dispersion, and have confirmed that the form of the con- 
densate wavefunction and density in the non-homogeneous phase is sensitive to the 
form of the repulsion. Work remains to be done regarding the effects of competing 
repulsive terms, in which analytic approaches may prove helpful. 
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